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Abstract. In this paper, we propose a general framework to design asymptotic preserving schemes 
for the Boltzmann kinetic kinetic and related equations. Numerically solving these equations are 
challenging due to the nonlinear stiff collision (source) terms induced by small mean free or relaxation 
time. We propose to penalize the nonlinear collision term by a BGK-type relaxation term, which 
can be solved explicitly even if discretized implicitly in time. Moreover, the BGK-type relaxation 
operator helps to drive the density distribution toward the local Maxwellian, thus natually imposes 
an asymptotic-preserving scheme in the Euler limit. The scheme so designed does not need any 
nonlinear iterative solver or the use of Wild Sum. It is uniformly stable in terms of the (possibly 
small) Knudsen number, and can capture the macroscopic fluid dynamic (Euler) limit even if the 
small scale determined by the Knudsen number is not numerically resolved. It is also consistent 
to the compressible Navicr-Stokes equations if the viscosity and heat conductivity are numerically 
resolved. The method is applicable to many other related problems, such as hyperbolic systems with 
stiff relaxation, and high order parabilic equations. 
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1. Introduction 

The Boltzmann equation describes the time evolution of the density distribution of a dilute gas of 
particles when the only interactions taken into account are binary elastic collisions. For space variable 
x E fi € R dx , particle velocity v E R d " (d v > 2), the Boltzmann eqaution reads: 

(1-1) |+^VJ=-Q(/) 

at e 

where / := f(t,x,v) is the time-dependent particles distribution function in the phase space. The 
parameter e > is the dimensionless Knudsen number which is the ratio the mean free path over a 
typical length scale such as the size of the spatial domain, thus measures the rarefiedness of the gas. 
The Boltzmann collision operator Q is a quadratic operator, 

(1.2) Q(f)(v)= f [ B(\v-v*\,cos0) (/*/'-/*/) dadv*. 
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We used the shorthanded notation / = f(v), /* = /(i>*), /' = f(v'), fi — /(«*)• The velocities of the 
colliding pairs an d (u',t>*) are related by 

■/ = V- -((« - u*) — |« - u*| cr), 

< =u - -((«-«*) + |«- «*k), 

with a S § The collision kernel 5 is a non- negative function which by physical arguments of 

invariance only depends on \v — u*| and cos# = u ■ a (where u = (u — — u*| is the normalized 

relative velocity). In this work we assume that B is locally integrable, given by 

B(|u|,cos0) = C 1 M 7 , 

for some 7 € (0, 1] and a constant C 7 > 0. 

Boltzmann's collision operator has the fundamental properties of conserving mass, momentum and 
energy: at the formal level 



(1.3) 



e(/)0(«)d« = o, 0(«) = i,«,h s 



and it satisfies the well-known Boltzmann's theorem 



G(/)log(/)d»>o. 



The functional — //log/ is the entropy of the solution. Boltzmann's iJ theorem implies that any 
equilibrium distribution function, i.e., any function which is a maximum of the entropy, has the form 
of a local Maxwellian distribution 

\u — v\ 2 " 



M p , U)T (v) 



exp 



2T 



where p, u, T are the density, macroscopic velocity and temperature of the gas, defined by 



(1.4) 
(1.5) 



f(v) dv 



T = 



1 



M.p,u,T(v), u= — 
P 

1 



u f(v) dv = — 



v A4p, u, T(v) dv 



\u — v\ f{v) dv = 



d v p 



\u — v\ A4p,u,T(v)dv 



Therefore, when the Knudsen number e > becomes very small, the macroscopic description, which 
describe the evolution of averaged quantities such as the density p, momentum pu and temperature 
T of the gas, by fluid dynamics equations, namely, the compressible Euler or Navier-Stokes equations, 
become adequate. Mor specifically, i.e. as e — > 0, the distribution function will converge to a local 
Maxwellian Ai, and the system (|1.2p becomes a closed system for the 2 + d v moments. The conserved 
quantities satisfy the classical Euler equations of gas dynamics for a mono-atomic gas: 

dp 



(1.6) 



dt 
dpu 

dE 
~dt 



X7 X ■ pu = 0, 
4- • {pu® u + pi) = 0, 



+ V x ■ {{E + p)u) 



0. 



where E represents the total energy 



E = -pu 2 
2 H 



— pT , 
2 y ' 



and I is the identity matrix. These equations constitute a system of 2 + d v equations in 3 + d v 
unknowns. The pressure is related to the internal energy by the constitutive relation for a polytropic 
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where the polytropic constant 7 = (d v + 2)/d v represents the ratio between specific heat at constant 
pressure and at constant volume, thus yielding p = pT. For small but non zero values of the Knudsen 
number e, the evolution equation for the moments can be derived by the so-called Chapman-Enskog 
expansion [8], applied to the Boltzmann equation. This approach gives the Navier-Stokes equations 
as a second order approximation with respect to e to the solution to the Boltzmann equation: 

dp e 



(1.7) 



dt 
dt 



x ■ Pe U s = 0, 
V s • {p e U £ (g)U £ + p e l) = eV x - [p £ <j(u £ )}, 



dE 

— j- + V x • (E £ + p £ ) u e ) = sV x ■ (p £ a(u £ ) u + k £ V X T £ ) . 



dt 

In these equations o~(u) denotes the strain-rate tensor given by 

a(u) = V x u + (V x u) T - — V x -ul 

while the viscosity p £ = p{T £ ) and the thermal conductivity k £ = k(T £ ) are defined according the 
linearized Boltzmann operator with respect to the local Maxwellian [TJ. 

The connection between kinetic and macroscopic fluid dynamics results from two properties of the 
collision operator: 

(i) conservation properties and an entropy relation that imply that the equilibria are Maxwellian 

distributions for the zeroth order limit ; 
(ii) the derivative of Q(f) satisfies a formal Fredholm alternative with a kernel related to the 

conservation properties of (i). 

Past progress on developing robust numerical schemes for kinetic equations that also work in the 
fluid regimes has been guided by the fluid dynamic limit, in the framework of asymptotic-preserving 
(AP) scheme. As summarized by Jin [31] . a scheme for the kinetic equation is AP if 

• it preserves the discrete analogy of the Chapman-Enskog expansion, namely, it is a suitable 
scheme for the kinetic equation, yet, when holding the mesh size and time step fixed and 
letting the Knudsen number go to zero, the scheme becomes a suitable scheme for the limiting 
Euler equations 

• implicit collision terms can be implemented explicitly, or at least more efficiently than using 
the Newton type solvers for nonlinear algebraic systems. 

To satisfy the first condition for AP, the scheme must be driven to the local Maxwellian when 
e — > 0. This can usually be achieved by a backward Euler or any L-stable ODE solvers for the 
collision term |32j . Such a scheme requires an implicit collision term to gaurantee a uniform stability 
in time. However, how to invert such an implicit, yet nonlocal and nonlinear, collision operator is a 
delicate numerical issue. Namely, it is hard to realize the second condition for AP schemes. 

Comparing with a multiphysics domain decomposition type method [H [132 El [29, 39, 43 ], the AP 
schemes avoid the coupling of physical equations of different scales where the coupling conditions are 
difficult to obtain, and interface locations hard to determine. The AP schemes are based on solving 
one equation- the kinetic equation, and they become robust macroscopic (fluid) solvers automatically 
when the Knudsen number goes to zero. An AP scheme implying a numerical convergence uniformly in 
the Knudsen number was proved by Golse-Jin-Levermore for linear transport equation in the diffusion 
regime 27 . This result can be extended to essentially all AP schemes, although the specific proof is 
problem dependent. For examples of AP schemes for kinetic equations in the fluid dynamic or diffusive 
regimes see for examples [111 [SJ [311 1311 13H I3H1 HH HI US] ■ The AP framework has also been extended 
in [T31 [TJ] for the study of the quasi-neutral limit of Euler-Poisson and Vlasov-Poisson systems, and 
in [TBI HO] for all-speed (Mach number) fluid equations bridging the passage from compressible flows 
to the incompressible flows. 

Since the Boltzmann collision term Q needs to be treated implicitly, how to invert it numerically 
becomes a tricky issue. One solution was offered by Gabetta, Pareschi and Toscani [25]. They first 
penalize Q by a linear function A/, and then absorb the linearly stiff part into the time variable to 
remove the stiffness. The remaining implicit nonlinear collision term is approximated by finite terms 
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in the Wild Sum, with the infinite sum replaced by the local Maxwellian. This yields a uniformly 
stable AP scheme for the collision term, capruring the Euler limit when e — > 0. Such a time-relaxed 
method was also used to develop AP Monte Carlo method, see [SHE]- Nevertheless, it seems that 
this method is not able to capture the compressible Navier-Stokes asymptotic for small e. 
When the collision operator Q is the BGK collision operator 



it is well-known that even an implicit collision term can be solved explicitly, using the property that 
Q preserves mass, momentum and energy. Our new idea is this paper is to utilize this property, and 
penalize the Boltzmann collision operator Q by the BGK operator: 



where A is the largest spectrum of the linearized collision operator of Q around the local Maxwellian M. . 
Now the first term on the right hand side of (II. 9p is either not stiff, or less stiff compared to the second 
term, thus can be discretized explicitly, so as to avoid inverting the nonlinear operator Q. The second 
term on the right hand side of (|1.9p is stiff, thus will be treated implicitly. Despite this, as mentioned 
earlier, the implicit BGK operator can be inverted explicitly. Therefore we arrive at a scheme which 
is uniformly stable in e, with an implicit source term that can be solved explicitly. In other words, in 
terms of handling the stiffness, the general Boltzmann collision operator can be handled as easily as 
the much simpler BGK operator, thus we significantly simplies an implicit Boltzmann solver! 

Although a linear penalty (by removing Ai on the right hand size from (|1.9p can also remove the 
stiffness, it does not have the AP property, unless one follows the Wild Sum procedure of [35]. The 
BGK operator that we use in (|I .9|) helps to drive / into M, thus preserves the Euler limit. This 
will be proved asymptotically for prepared initial data (namely data near Ai), and demonstrated 
numerically even for general initial data. Moreover, we will prove asymptotically that, for suitably 
small time-step, this method is also consistent to the Navier-Stokes equations q 1.7(1 for e << I. 

Our method is partly motivated by the work of Haack, Jin and Liu [3D] , where by subtracting the 
leading linear part of the pressure in the compresible Euler equations with a low Mach number, the 
nonlinear stiffness in the pressure term due to the low Mach number is removed and an AP scheme 
was proposed for the compressible Euler or Navier-Stokes equations that capture the incompressible 
Euler or Navier-Stokes limit when the Mach number goes to zero. 

Our method is not restricted to the Boltzmann equation. It applies to general nonlinear hyperbolic 
systems with stiff nonlinear relaxation terms [Till GH3 H2J [TJJ , and higher-order parabolic equations 
(see section 5). Moreover, it applies to any stiff source term that admits a stable local equilibrium. 

In the following sections, we present a class of asymptotic preserving schemes designed for kinetic 
equations even if the general framework can be applied to other partial differential equations. We 
will focus on the Boltzmann equation and its hydrodynamic limit. We present different numerical 
tests to illustrate the efficiency of the present method. We treat particularly a multi-scale problem 
where the Knudsen numer e depends on the space variable and takes different values ranging from 
(hydrodynamic regime) to one (kinetic regime). Finally, the last part is devoted to the design of 
numerical schemes for nonlinear Fokker-Planck equations for which the asymptotic preserving scheme 
can be used to remove the CFL constraint of a parabolic equation. 



Since out method does not depend on the discretization of the spatial derivative, but only on the 
structure of the stiff source term, we will first present in in the simplest framework for stiff ordinary 
differential equations. 

Let us consider a Hilbert space H and the following nonlinear autonomous ordinary differential 

system 



(1.8) 



Qbgk = M-f, 



(1.9) 



Q = [Q - X(M - /)] + X[M - f] 



2. An Asymptotic Preserving (AP) stiff ODE solver 
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• the solution to (|2.1[) converges to the steady state M. when time goes to infinity, and the 
spectrum of VQ(/) C C~ = z£ Cr,Im(z) < 0, 

< a < ||Ve(/)|| < £, V/ € ff. 

Remark 2.1. The second hypothesis above is certainly not the most general, but is convenient for our 
purpose. The lower bound implies that the solution converges to the steady state M, while the upper 
bound is a sufficient condition for existence and uniqueness of a global solution. 

When e becomes small, the differential equation (|2.1|) becomes stiff and explicit schemes are subject 
to severe stability constraints. Of course, implicit schemes allow larger time step, but new difficulty 
arises in seeking the numerical solution of a fully nonlinear problem at each time step. Here we want 
to combine both advantages of implicit and explicit schemes : large time step for stiff problems and 
low computational cost of the numerical solution at each time step. 

We denote by /" an approximation of f(t n ) with t™ = n At and the time step At > 0, Two classical 
procedures handle the aforementioned difficulties well. One is to linearize the unknown Q(f n+1 ) at 
time step t n+1 around / at the previous time step f n : 

Q(f n+1 ) « Q(P) + VQ(/")(/ n+1 - f n ) 

yielding a problem that only needs to solve a linear system with coefficient matrix depending on 
VQ(/ n ) [H|. This approach gives a uniformly stable time discretization without nonlinear solvers, 
however, it is not AP since the right hand size, as e — * 0, does not project f n+1 to the local equilibrium 
/ = M., even if f n = M. . The second approach, introduced by [2Sj , takes 

Q(/) = [Q(/)-m/]+m/- 

As mentioned in the introduction, it uses the Wild Sum expansion for Q(f) on the right side, which 
is truncated and the remaining infinite series is replaced by the local Maxwellian in order to be AP 
for the Euler limit. 

Under our hypothesis, the asymptotic behavior of the exact solution f £ is known when e — > 0. 
Therefore, we split the source term of (|2.ip in a stiff and non- (or less) stiff part as 

Q(f) _ go) - pjf) , no 

see 



non stiff part stiff part 

where P(f) is a well balanced, i.e. preserving the steady state, P(A4) — 0, linear operator and is close 
to the source term Q(f). For instance, performing a simple Taylor expansion, we get 

Q(f) = Q(M) + VQ(M)(f-M) + 0(\\f-M\\ 2 H ) 

and we may choose 

P(f) := VQ(M)(f-M). 
Since it is not always possible to compute exactly VQ(A4), we may simply choose 

P(f) := L(f-M), 

where L is an estimate of \7Q(A4). 

Now, we simply apply a first order implicit-explicit (IMEX) scheme for the time discretization of 
(EU): 

- r _ qct) - p{p) p(f n+i ) 

[Z - Z) At e + e ' 

or 

f n+1 = [si - AtVQiM)]' 1 [sf n + At(Q(f n ) - P(P)) - AtVQ{M)M]. 

This method is easy to implement, since f n+1 is linear in the right hand side of l|2.2j) . For linear 
problems, we have the following result: 

Theorem 2.2. Consider the differential system A2. 1)) with Q(f) = —A/, where Re(X) > 0. Set 
P(f) := —v\f with v > 0. Then, the scheme 12. 8\) is A-stable and L-stable for v > 1/2. 
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Proof. For linear systems, the scheme simple reads 

+1 = e + (v — 1) A At / AAt \ 

7 e + vXAt J \ e + uXAt) 1 ' 

Observe that v = gives the explicit Euler scheme, which is stable only for At < e/X, whereas for 
< v < 1, it yields the so-called ^-scheme, which is A-stable for v > 1/2. For v = 1 it corresponds to 
the A-stable implicit Euler scheme. Moreover, for v > 1, the scheme is A-stable, that is 

XAt 
f ^ AAf 

where |1 — ~| < 1 for v > 1/2. This is also the condition for the L-stability [26]. □ 



iir +1 ii* < 



i - 



link- 




H) 







\f n \\ H for e~0, 



To improve the numerical accuracy, second order schemes are sometimes more desirable. Thus, 
we propose the following second order IMEX extension. Assume that an approximate solution /" is 
known at time t n , we compute a first approximation at time f n+1 / 2 = t n + At/2 using a first order 
IMEX scheme and next apply the trapezoidal rule and the mid-point formula. The scheme reads 

f* - r Q(f n ) - p(n , p(n 



(2.3) 



2- 

At e e 

/»+!-/» Q(P) - P(P) , P(f n ) + P(f n+1 ) 





At e 2 s 

Note that both (|2.2p and (|2.3p are AP for prepared initial data. Let us use (|2.3j) as the example. 
Assume f n = M + 0(e). Then Q{f n ) = 0{e),P(f n ) = 0(e), thus the first step in j23|) gives 
P(f*) = 0(e). This further implies that Q(f*) = 0(e). Applying all these in the second step of (|2.3p 
gives P(f n+1 ) = 0(e), thus f n+1 = M + 0(e), which is the desired AP property. 

To illustrate the efficiency of (I2.2[) and (|2.3[) in various situations, we consider a simple linear 
problem with different scales for which only some components rapidly converge to a steady state 
whereas the remaining part oscillates. We solve 

(2-4) Q(f) =Af, 

where 

(2.5) A = 

for which the eigenvalues are Sp(A) = { — 1000 + i, —1000 — i, i}. The first block represents the 
fast scales whereas the last one is the oscillating part. Indeed, the first components go to zero 
exponentially fast whereas the third one oscillates with respect to time with a period of 27r. We want 
to solve accurately the oscillating part with a large time step without resolving small scales. Then, 
we apply the first order (12. 2|) and second order (|2.3p schemes by choosing 

P(f) = vAf, 

with v > 0. Here we take a large time step At = 0.3 and v = 2, which means that P(f) has the 
same structure of Q(f) but the eigenvalues are over estimated. Thus, fast scales are under-resolved 
whereas this times step is a good discretization of the third oscillating component. Therefore, an 
efficient AP scheme would give an accurate behavior of the slow oscillating scale with large time step 
with respect to the fast scale. It clearly appears in Figure [T] that the time step is too large to give 
accurate results for the first order scheme (|2.2p : the solution is stable but the oscillation of the third 
component is damped for this time step which is too large. This approximation is compared with the 
one obtained with a first order explicit Euler using a times step ten times smaller. We also compare 
the numerical solution of the second order scheme (|2.3p with the one obtained using a second order 
explicit Runge-Kutta scheme corresponding to v = with a time step three hundred times smaller. 
In Figure [TJ we observe the stability and good accuracy of the second order scheme (|2.3|) . Note that 
for the same time step, the explicit Runge-Kutta scheme blows-up! 

In the following sections we apply this approach to the Boltzmann equation and verify its accuracy 
and efficiency on several classical problems dealing with fluid, kinetic and multi-scale regimes. 
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3. Application to the Boltzmann equation 

We now extend the stiff ODE solver of the previous section to the Boltzmann eqaution To 
this aim, we rewrite the Boltzmann equation (jl.ip in the following form 



(3.1) 



at e e 



where the operator P is a "well balanced relaxation approximation" of Q(f), which means that it 
satisfies the following (balance law) 

P(f)<f>(v)dv = 0, cf>(v) = l,v,\v\ 2 , 

and preserves the steady state i.e. P(A4 p , u ,t) — where M. p . u .t is the Maxwellian distribution 
associated to /?, u and T given by (jl.4|l . Moreover, it is a relaxation operator in velocity 

(3.2) P(f)=p[M p , u , T (v)-f(v)}. 

For instance, P(f) can be computed from an expansion of the Boltzmann operator with respect to 
Mp,u,T: 

Q(f) ~ Q{M p , u ,t) + VQ(M,, w ,r) [M P , U , T - f] ■ 

Thus, we choose (3 > as an upper bound of the operator VQ(M p , u ,t)- Then P(f) given by (|3.2|) is 
just the BGK collisional operator [3]. 

Since the convection term in (|3.1|) is not stiff, we will treat it explicitly. For source terms on the 
right hand side of (|3.ip will be handled using the ODE solver in the previous section. For example, if 
the first order scheme (12.21) is used, then we have 



(3.3) 



f n+l _ f n Q{r) _ p {fn) p^n+l) 

+ v ■ VxJ = 1 • 



At 

f°(x,v) = fo(x,v) . 
Using the relaxation structure of P(f), it can be written as 



r +1 - — : — 77T7 [/" - At » V,/"] + At ■ 



e + /5At J ' £ + (3At 

M r 



where is the Maxwellian distribution computed from f n+1 . 
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Although (|3.4p appears nonlinear ly implicit, it can be computed explicitly. Specifically, upon 
multiplying (|3.4p by (j){v) defined in ()1.3|) . and use the conservation property of Q and P and the 
definition of M in l|1.4j) . one gets 



U 



n+1 



e + f3At 



or simply 



U 



n+1 



<t>(r-Atvv x f n )dv+ 



4(r-Atvv x p)dv. 



(3 At 



U 



n+1 



Thus U n+1 can be obtained explicitly, which defines A4 n+1 . Now f n+1 can be obtained from 
explicitly. In summary, although (|3.3p is nonlinearly implicit, it can be solved explicitly, thus satisfies 
the second condition of an AP scheme. 

We define the macroscopic quantity U by U :— (p,pu,T) computed from /. Clearly, the scheme 
satisfies the following properties 



Proposition 3.1. Consider the numerical solution given by i3.S\) . Then, 

(i) If e — > and /" = A4 n + 0(e), then the scheme &3. 3]) is asymptotic preserving, that is 
f n+1 = A4' n+1 + 0(e) ; thus the scheme is AP to the Euler limit in the sence that, when e — > 0, 
the (moments of the) scheme becomes a consistent discretization of the Euler system \1.6\) . 
(ii) For e < 1 and there exists a constant C > such that 



(3.4) 



jn-\-\ £fi 


+ 


Tjn+1 _ jjn 


At 




At 



<c, 



then the scheme k3.3\) asymptotically becomes a first order in time approximation of the com- 
pressible Navier- Stokes 

Proof. We easily first check that for e — > and /" = M n , we get f n+1 = A4 n+1 . Therefore, we 
multiply (|3.3p by (1, v, \v\ 2 /2) and integrate with respect to v, which yields that U n is given by a time 
explicit scheme of the Euler system (|1.6j) . 

Now let us prove (ii). We apply the classical Chapman-Enskog expansion: 

(3.5) /" = M n + eg n 

and integrate ()3.3j) with respect to v £ M. d " . By using the conservation properties of the Boltzmann 
operator (|1.3|) and of the well-balanced approximation P(f), 



(3.6) 



jjn+l _ jjn 
At 



V„ 



1 

V 

2 



v (M n + eg n )dv 



For eg — 0, this is the compressible Euler equations (|1.6p . Thus, a consistent approximation of 
the compressible Navier-Stokes is directly related to a consistent approximation of g n . Inserting 
decomposition (|3.5p into the scheme (|3.3p gives 

M n+l_ M n /f ,n+l 

- + vV x M n + e 



At 



At 



W x g n 



Q(M n + eg r - 



- [i3( P n )g n - P( P n+1 )g n+1 ] 



Since Q is a bilinear and Q(A4) = 0, one has 

Q(M+eg) = Q(M) + e £ M (g) + e 2 Q(g), 
where Cm is the linearized collision operator with respect to M.. Thus, we get 

M n+1 - M" 



+ 



At 

n+1 _ 



a" 



(3.7) 



At 

C M {g n ) - vV x M n 



- W l )a n - P(p n+1 )g n+1 ] 

vV x g n -Q{g n ) 
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It is well known that Cm is a non-positive self-adjoint operator on L M denned by the set 



{<p: <pM- 1/2 e L 2 (K d ")} 



and that its kernel is N{Cm) = Span{.M, v A4, \v\ 2 A4}. Let Hm be the orthogonal projection in L 
onto M{Cm)- After easy computations in the orthogonal basis, one finds that 



M 



where 



n^W = — 
p 



niQ = I ip dv, mi 



too 



v — u ( \v — u\ 2 d l 



T 



-777-1 



2T 2 1 '"' 



2T 2 ' 



(v — u)ip dv, TO 2 
It is easy to verify that ILm" (M n ) = M n and 

iWg") = n M n(g n+1 ) = u M n(Q(g n )) = (<?")) = o 

Then applying the orthogonal projection I — Um" to (|3.7p . it yields 

' M n+1 - M n 



(i-n 



g n+l _ g n 



At 



+ (/3(p n )g n - (3( P n+1 )g n+1 ) 



At 



(l-U M n)(vV x g n )-Q(9 n ) 



= C M (g n )~(l-U Mn )(vV x M n ). 
Finally, it remains to estimate 

(I - H M n 

Using a Taylor expansion we find that 

M n+1 = M n 



M n+1 - M n 
At 



1 P n+1 -P n V- u n n+1 n> 
1 + — + — (u + - u n ) 



v-u n \ 2 d\ T n+1 — T n 



2T n 



0{At 2 ) 



and by definition of 11^ 

n M n(M n+1 ) 



M" i 1 + pn+1 pn + V -^~ (u n+1 - u n ^ 



T" 



M' 



n (\ v~u n \ 2 _ d 
2T n 2 



\v-u n \ 2 d\ T n+1 -T 
2T n 

n+l 



T" 



p n T n (P n+1 ~ + ^r("" +1 " «")" 



v — ii n n n+1 — n n 

+ M» K +1 - «") + 0(At 2 ). 

Thus, under the assumption (|3.4p , we have 



(i-iW.) 



(M n+1 - M n 



V At 

and the residual distribution function is given by 



0(At) 



g n = ( (I - U M n) (v ■ \7 x M n ) ) + 0(e) + 0(At). 



Now, substituting this latter expression in (|3.6[) . we get 

/ v 



U n+1 - U n 
At 



V X -F(U) = -eV x 

+0(eAt + e 2 



V Q9 V 
U,|2 



C^ n ((ld-H M «)(vV x M n ))dv 



\ 



where 



pu 

F(U) = | pwg/u+pl 
(E + p)u 
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To complete the proof, it remains to compute the term in 0(e). An easy computation first gives 



(i-rwO(«. v x M n ) = 



B 



Vu+ (Vu) T - -V- ul 



A 



VT 



M(v), 



with 



.4 



d + 2 



2T 



T 



B =2 



1 / (v — u) ® (v — u) \i 



2T 



Therefore, it yields 

Cjj n ((I-IL M n)(vV x M n )) 



^»(BAf) ( Vu+ (Vu) T 
VT 



dT 



-V-uI 



VT' 



Substituting this expression in (|3.6| . we get a consistent time discretization scheme to the compressible 
Navier-Stokes system where the term of order of e is given by 





Me 0"(«e) 

H e cr{u e ) u 



with 



a(u) = V x u + (y x u) J 



2 



rf,, 



V, 



' I 



while the viscosity fi £ = /u(T e ) and the thermal conductivity k £ = k(T £ ) are defined according the 
linearized Boltzmann operator with respect to the local Maxwellian pQ. □ 

Remark 3.2. To capture the Navier-Stokes approximation that has 0(e) viscosity and heat conduc- 
tivity, one needs the mesh size and cAt to be o(e) (c is a characteristic speed). Thus conclusion (ii) 
in the above proposition shows that the scheme is consistent to the Navier-Stokes equations provided 
that the viscous terms are resolved, while to capturing the Euler limit one can use mesh size and cAt 
much larger than e, in the usual sense of asymptotic-preserving. 

4. Numerical tests 

In this section we perform several numerical simulations for the Boltzmann equation in different 
asymptotic regimes in order to check the performance (in stability and accuracy) of our methods. We 
have implemented the first order (|2.2p and second order ()2.3|) scheme for the approximation of the 
Boltzmann equation. Here, the Boltzmann collision operator is discretized by a deterministic method 
[HI [19l [2Qj [22] , which gives a spectrally accurate approximation. A classical second order finite volume 
scheme with slope limiters is applied for the transport operator. 

4.1. Approximation of smooth solutions. This test is used to evaluate the order of accuracy of 
our new methods. More precisely, we want to show that our methods (|2.2p and (|2.3p are uniformly 
accurate with respect to the parameter e > 0. We consider the Boltzmann equation (11. ip in 14x24- 
We take a smooth initial data 



fo(x,v) 



Po(x) 



exp 



(x,v) € [—L,L] x 



2nT (x) "V 2T (x), 

with po(x) = (11 — 9tanh(x))/10, T (x) = (3 — tanh(x))/4, L — 1 and assume specular reflection 
boundary conditions in x. Numerical solutions are computed from different phase space meshes : the 
number of point in space is n x — 50, 100, 200,. ..,1600 and the number of points in velocity is with 
n v = 8,.. .,64 (for which the spectral accuracy is achieved), the time step is computed such that the 
CFL condition for the transport is satisfied At < Ax/v max , where Aa; is the space step and u max = 7 
is the truncation of the velocity domain. Then different values of e are considered starting from the 
fully kinetic regime e — 1, up to the fluid limit e = 10~ 5 corresponding to the solution of the Euler 
system (11.61) . The final time is T max = 1 such that the solution is smooth for the different regimes. 
An estimation of the relative error in L p norm is given by 

\fh(t) ~ f2h(t)\\ P y 



£2h 



max 

te(o,T) 



\h\\ P 



1 < P < +oo, 
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where fh represents the approximation computed from a grid of order h. The numerical scheme is 
said to be fc-th order if e-ih < C h k , for all < h <C 1. 



£= 1.E-5 

e= 1.E-2 
e= 1.E-1 

£ = 2.E-1 
I E=5.E-1 

e= 1.E+0 



-2 




1.9 x 



log(h) 

(1) 



-0.6 



-0.4 



- e = 1.E-5 


— I 






e = 1.E-2 








e = 1 .E-1 


* 






e = 2. E-1 


— B 






e =5.E-1 


------ 






8 = 1.E+0 




















1.9 x 






-'W 












□ 









-1.6 -1.4 -1.2 -1 -0.! 

log(h) 

(2) 



-0.6 



-0.4 



Figure 2. The L 1 and L°° errors of the second order method l|2.3p for different 
values of the Knudsen number s = 1CP 5 ,...,1. 

In Figure [2j the L 1 and L°° errors of the second order method (|2.3|) are presented. They show a 
uniformly second order convergence rate (an estimation of the slope is 1.9) in space and time (the 
velocity discretization is spectrally accuracy in v thus does not contribute much to the errors). The 
time step is not constrained by the value of e, showing a uniform stability in time. 

4.2. The Sod tube problem. This test deals with the numerical solution to the ld x x 2d v Boltzmann 
equation for Maxwellian molecules (7 = 0). We present numerical simulations for one dimensional 
Riemann problem and compute an approximation for different Knudsen numbers, from rarefied regime 
to the fluid regime. 

Here, the initial data corresponding to the Boltzmann equations are given by the Maxwellian 
distributions computed from the following macroscopic quantities 

{puUuTfi = (1,0,1), if0<x<0.5, 

(p r , u r , T r ) = (0.125, 0, 0.25) , if 0.5 < x < 1 . 

We perform several computations for e = 1, 10 _1 , 10~ 2 ,...,10~ 4 . In Figure 02 we only show the 
results obtained in the kinetic regime (10 -2 ) using a spectral scheme for the discretization of the 
collision operator |22j (with n v = 32 2 and a truncation of the velocity domain u max = 7) and second 
order explicit Runge-Kutta and second order method (|2.3p for the time discretization with a time step 
At = 0.005 satisfying the CFL condition for the transport part (with n x — 100). For such a value of 
£, the problem is not stiff and this test is only performed to compare the accuracy of our second order 
scheme (|2.3p with the classical (second order) Runge-Kutta method. We present several snapshots of 
the density, mean velocity, temperature and heat flux 

Q(t,x) := - / (v - u e ) \v - u e \ 2 f e (t,x,v)dv 

at different time t = 0.10 and 0.20. Both results agree well with only n x — 100 in the space domain 
and n v — 32 for the velocity space. Thus, in the kinetic regime our second order method (|2.3p gives 
the same accuracy as a second order fully explicit scheme without any additional computational effort. 

Now, we investigate the cases of small values of e for which an explicit scheme requires the time 
step to be of order O(e). In order to evaluate the accuracy of our method (|2.3[) in the Navier-Stokcs 
regime (for small e -C 1 but not negligible), we compared the numerical solution for e = 10~ 3 with 
one obtained with a small time step At = 0(e) (for which the computation is still feasible). Note 
that a direct comparison with the numerical solution to the compressible Navier-Stokes system f|l ,7|) 
is difficult since the viscosity fi s — fi(T e ) and the thermal conductivity k £ — k(T s ) are not explicitly 
known. Therefore, in Figure|H we report the numerical results for e = 10~ 3 and propose a comparison 
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Figure 3. Sod tube problem (e = 1CP 2 ), dots (x) represent the numerical solution 
obtained with our second order method (|2.3|) and lines with the Runge-Kutta method: 
evolution of (1) the density p, (2) mean velocity u, (3) temperature T and (4) heat 
flux Q at time t = 0.05, 0.1, 0.15 and 0.2. 




between the numerical solution obtained with the scheme (|2.3j) and the one obtained with a second 
order explicit Runge-Kutta method. In this case, the behavior of macroscopic quantities (density, 
mean velocity, temperature and heat flux) agree very well even it the time step is at least ten times 
larger with our method (|2.2|) or (|2.3|) . 

Finally in Figure [HI we compare the numerical solution of the Boltzmann equation (jl.ip with the 
numerical solution to the compressible Navier-Stokes system derived from the BGK model since the 
viscosity and heat conductivity are in that case explicitly known [2] . To approximate the compressible 
Navier-Stokes system, we apply a second order Lax-Friedrich scheme using a large number of points 
(n x = 1000) whereas we only used n x = 100, and 200 points in space and n 2 , = 32 2 points in velocity 
for the approximation of the kinetic equation (jl.lj) . In this problem, the density, mean velocity and 
temperature are relatively close to the one obtained with the approximation of the Navier-Stokes 
system. Even the qualitative behavior of the heat flux agrees well with the heat flux corresponding to 
the compressible Navier-Stokes system k e V x T e , with k e — p e T £ (see Figured]), yet some differences 
can be observed, which means that the use of BGK models to derive macroscopic models has a strong 
influence on the heat flux. 



4.3. A problem with mixing regimes. Now we consider the Boltzmann equation (|1.1|) with the 
Knudsen number e > depending on the space variable in a wide range of mixing scales. 



ASYMPTOTIC PRESERVING SCHEMES FOR KINETIC EQUATIONS 



13 




(1) 





0.2 



(2) 




(3) 



(4) 



Figure 4. Sod tube problem (e = 1CP 3 ), dots (x) represent the numerical solution 
obtained with our second order method (|2.3|) and lines with the Runge-Kutta method: 
evolution of (1) the density p, (2) mean velocity u, (3) temperature T and (4) heat 
flux Q at time t = 0.05, 0.1, 0.15 and 0.2. 



This kind of problem was already studied by several 
authors for the BGK model [T7] or the radiative transfer 
equation [35]. In this problem, £ : Rh M + is given by 

e(x) = £ + - [tanh(l -Ux) + tanh(l + 11 x) ] , 
which varies smoothly from Eq to 0(1). 



This numerical test is difficult because different scales are involved. It requires a good accuracy 
of the numerical scheme for all range of s. In order to focus on the multi-scale nature we only 
consider periodic boundary conditions, even if the method has also been used with specular reflection 
in space. Furthermore, to increase the difficulty we consider an initial data which is far from the local 
equilibrium of the collision operator: 

x G [~L,L], v G M 2 
with u = (3/4,-3/4), 

. , 2 + sin(fc:c) 5 + 2cos(fca;) 
Po( x ) = g ' °( X ) = 20 
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Figure 5. Sod tube problem (e = 1CP 4 ), comparison of the numerical solution to the 
Boltzmann equation with the our second order method (|2.3|) represented with dots 
(x) with the numerical solution to the compressible Navier-Stokes system (lines): 
evolution of (1) the density p, (2) mean velocity u, (3) temperature T and (4) heat 
flux Q at time t = 0.05, 0.1, 0.15 and 0.2. 



where k — ir/L and L = 1/2. 

Here we cannot compare the numerical solution with the one obtained by a macroscopic model. 
From the numerical simulations, we observe that the solution is smooth during a short time and some 
discontinuities are formed in the region where the Knudsen number e is very small and then propagate 
into the physical domain. 

On the one hand, we only take £o = 10 -3 in order to propose a comparison of numerical solutions 
computed with a second order method using a time step At = 0.001 (such that the CFL condition 
for the transport part is satisfied) and the one by the second order explicit Runge-Kutta method 
with a smaller time step At — 0.0001) to get stability. The number of points in space is n x = 200 
and in velocity is nf, = 32 2 . Clearly, in Figure the results are in good agreement even if our new 
method does not solve accurately small time scales when the solution is for from the local equilibrium. 
Moreover in Figure [71 we present numerical results with only n x = 50 and n x — 200, and v? v = 32 2 
to show the performance of the method with a small number of discretization points in space. With 
n x — 50 points the qualitative behavior of the macroscopic quantities (p, u, T) is fairly good. 

On the other hand, we have performed different numerical results when Eq = 10~ 4 , then the 
variations of e starts from 10~ 4 to 1 in the space domain. In that case, the computational time of 
a fully explicit scheme would be more than one hundred times larger than the one required for the 
asymptotic preserving scheme (|2.3|) . We observe that discontinuities appear on the density, mean 
velocity and temperature and then propagate accurately into the domain. The shock speed is roughly 
the same for the different numerical resolutions. Therefore, this method gives a very good compromise 
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between accuracy and stability for the different regimes. Numerical results are not plotted since they 
are relatively close to the ones presented in Figures E] and [71 
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Figure 6. Mixing regime problem (eq = 10~ 3 ), comparison of the numerical solution 
to the Boltzmann equation with the second order method (|2.3|) represented with dots 
(x) with the numerical solution obtained with the explicit Runge-Kutta method using 
a small time step (line): evolution of (1) the density p, (2) mean velocity u, (3) 
temperature T at time t — 0.25, 0.5 and 0.75. 



5. Other applications: numerical stability 

In this section, we want to illustate the efficiency of the asymptotic preserving scheme to treat high 
order differential operators. We have already applied such a scheme for Willmore flow (fourth order 
differential operator Here, we consider the flow of gas in a two dimensional porous medium 

with initial density go(v) > 0. The distribution function g(t, v) then satisfies the nonlinear degenerate 
parabolic equation 

(5.1) §=A^, 
where m > 1 is a physical constant. Assuming that 

/ (I + \v\ 2 ) g {v)dv < +oo, 
Jr 2 

J. A. Carrillo and G. Toscani [7 proved that g(t,v) behaves asymptotically in a self-similar way like 
the Barenblatt-Pattle solution, as t — » +oo. More precisely, it is easy to see that if we consider the 
change of variables 

(5.2) ^^^K l0g(S(<)) '^))' 
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Figure 7. Mixing regime problem (eq = 10~ 3 ), comparison of the numerical solution 
to the Boltzmann equation obtained with the AP scheme (|2.3|) using n x = 50 (dots 
x) and n x — 200 points (line): evolution of (1) the density p, (2) mean velocity u, (3) 
temperature T at time t = 0.25, 0.5 and 0.75. 



where s(t) := \J\ + 2i, the new distribution function / is solution to 

^ = V„>/ + V„/ m ), 
and converges to the Barenblatt-Pattle distribution 

m — 1 , 



M{v) 



C 



2 m 



i/(m-i) 



where C is uniquely determined and depends on the initial mass go but not on the "details" of the 
initial data. 

Instead of working on (|5.ip directly, we will study the asymptotic decay towards its equilibrium. 
The key argument on the proof of J. A. Carrillo and G. Toscani is the control of the entropy functional 



which satisfies 



H(f) 

dH(f) 
dt 



\v\ 2 f(t,v) 



f(t,v) 



-f m (t,v) 



m — 1 



dv. 



dv<0 



or the control of the relative entropy H(f\A4) = H(f) — H(M) with respect to the steady state Ai. 

Numerical discretization of this problem leads to the following difficulty : explicit schemes are 
constrained by a CFL condition At ~ Av 2 whereas implicit schemes require the numerical resolution 
of a nonlinear problem at each time step (with a local constraint on the time step). We refer to [51I2T] 
for a fully implicit approximation preserving steady states for nonlinear Fokker-Planck type equations. 
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Here we do not focus on the velocity discretization, but only want to apply our splitting operator 
technique to remove this severe constraint on the time step. Here the parameter e does not represent 
a physical time scale but is only related to the velocity space discretization Av. Therefore, we set 
Q(f) = V„ • (v f + V„/ m ) and P(f) = VQ(M) (f-M), which leads to the following decomposition 

% = A v (r - mM m - l f) + V v ■ (vf + mV„ (M" 1 ' 1 /)) . 

Ot v v > s „ < 

non stiff part stiff linear part 

Then we apply a simple IMEX scheme which only requires the numerical resolution of a linear system 
at each time step. 

We choose m = 3 and a discontinuous initial datum far from the equilibrium 

fo(v)= Yj Jfi 1 B(0,ro)( v ~ v k,l) 

(e{i,2} fee{o,...,n-i} 

where n — 12, tq = 1/4 and Vk,i — le l9k , with 9k = 2kir/n, k — 0, ... ,n — 1. We use a standard 
velocity discretization in the velocity space based on an upwind finite volume approximation for the 
transport term and a center difference for the diffusive part. We take n 2 — 120 2 in velocity and a 
time step At = 0.02 which is much larger than the time step satisfying a classical CFL condition for 
this problem At ~ 0(Av 2 ). The numerical scheme (|2.2|) is still stable and the numerical solution 
preserves nonegativity at each time step (see Figure [5])! For large time, the solution converges to an 
approximation of the steady state even if the present scheme is not exactly well-balanced (it does not 
preserve exactly the steady state). Moreover, to get a better idea on the behavior of the numerical 
solution, we plot the evolution of the entropy and its dissipation for different time steps. More 
surprisingly, the numerical entropy is decreasing and the dissipation converges towards zero when 
times goes to infinity. 



6. Conclusion 

We have proposed a new class of numerical schemes for physical problems with multiple time and 
spatial scales described by a still nonlinear source term. A prototype equation of this type is the 
Boltzmann equation for rarified gas. When the Knudsen number is small, the stiff collision term of 
the Boltzmann equation drives the density distribution to the local Maxwellian, thus the macroscopic 
quantities such as mass, velocity and temperature are be evolved according to fluid dynamic equations 
such as the Euler or Navier-Stokes equations. Asmptotic-preserving (AP) schemes for kinetic equations 
have been successful since they capture the fluid dynamic behavior even without numerically resolving 
the small Knudsen number. However, the AP schemes need to treat the stiff collision terms implicitly, 
thus it yeilds a complicated numerical algebraic problem due to the nonlinearity and nonlocality of 
the collision term. In this paper, we propose to augement the nonlinear Boltzmann collision operator 
by a much simpler BGK collision operator, and impose implicity only on the BGK operators which 
can be handled much more easily. We show that this method is AP in the Euler regime, and is also 
consistent to the Navier-Stokes approximations for suitably small time steps and mesh sizes. Numerical 
examples, including those with mixing scales and non-local-Maxwellian initial data, decomstrate the 
AP property as well as uniform convergence (in the Knudsen number) of this method. 

This method can be extended to a wide class of PDEs (or ODEs) with stiff source terms that admit 
a stable and unique local equilibirum. We use the Fokker-Planck equation as an example to illustrate 
this point, and will pursue more applications in the future. 



Acknowledgments. F. Filbet thanks Ph. Laurengot, M. Lemou, P. Degond and L. Pareschi for 
usefull discussions on the topic. 



18 



FRANCIS FILBET AND SHI JIN 



level set of numerical solution, t = 0.1 level set of numerical solution, t = 0.4 




Figure 8. Nonlinear Fokker-Planck solution: convergence toward equilibrium 
(Barenblatt-Pattle distribution) obtained with the first order method (|2.2p using 
n x = 100 at time t — 0.1, 0.4, 0.8, 1.0, 1.2 and 4 with a large time step. 
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Figure 9. Nonlinear Fokker-Planck solution: convergence toward equilibrium 
(Barenblatt-Pattle distribution) obtained with the first order method (|2.2[) using 
n x = 100 with At = 0.02 and 0.001. 
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